home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / LIBIP / SGL_STAT.C < prev    next >
C/C++ Source or Header  |  1999-09-11  |  8KB  |  303 lines
  1. /* 
  2.  * sgl_stat.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * SGL_STATS
  11.  *
  12.  * determine geometrical properties of a cluster of linear segments
  13.  * in a given SGL
  14.  */
  15.  
  16. #include <stdio.h>
  17. #include <string.h>
  18. #include <malloc.h>
  19. #include <math.h>
  20. #include "ip.h"
  21.  
  22. #define    TWO        2
  23. #define    THREE        3
  24. #define    DIM         TWO
  25.  
  26. #define    SN_LOLIM    0.002          /* lower limit: 1/512 */
  27. #define    CS_LOLIM    0.002
  28.  
  29. #undef DEBUG
  30.  
  31. /*
  32.  * nematic_op()
  33.  *   DESCRIPTION:
  34.  *     nematic_op evaluates nematic order parameter S = <3*cos_th*cos_th - 1>/2 (3d)
  35.  *     or S = <2*cos_th*cos_th - 1> (2d) for the ensemble of clusters (SGLs);
  36.  *   ARGUMENTS:
  37.  *     n_sgl(int) number of segment group lists
  38.  *     n_segm(int) number of clusters
  39.  *     dirns(float *) cluster orientations
  40.  *   RETURN VALUE:
  41.  *     total angular bend for the contour
  42.  *   COMMENTS:
  43.  *     note: *dirns is a ptr to an array containing the slopes of the
  44.  **             linear segments in a given SGL
  45.  **        loop over the n_sgl <= n_segm of non_empty SGLs only:
  46.  **             *dirns is expected to be initialized to 999.0; this
  47.  **             serves to identify unused elements; 
  48.  **             (see also: slope() in xsgll.c);
  49.  **        for each SGL (cluster) the array *dirns holds the average slope 
  50.  **             of the segments in that cluster; since there is no preferred
  51.  **             orientation of segments, corresponding directions are
  52.  **             defined to within modulo PI; it is with respect to the
  53.  **             the average direction that the order parameter S is
  54.  **             evaluated 
  55.  */
  56.  
  57. float
  58. nematic_op (int n_sgl, int n_segm, float *dirns)
  59. {
  60.   float slpi, av_slope;
  61.   double sn_2th, cs_2th, th;
  62.   double cs2_th, tg_th, s;
  63.   float order_parm;
  64.   int i;
  65.  
  66. #ifdef DEBUG
  67.   printf ("\n...nematic order parameter for %d clusters:\n", n_segm);
  68.   printf ("...cluster orientations:\n");
  69.   for (i = 0; i < n_segm; i++) {
  70.     if (*(dirns + i) < 999.0)
  71.       printf ("    SGL<%d>: %6.2f\n", i, *(dirns + i));
  72.   }
  73. #endif
  74.  
  75. /* determine preferred direction */
  76.   sn_2th = cs_2th = 0.0;
  77.   for (i = 0; i < n_segm; i++) {
  78.     if ((slpi = *(dirns + i)) < 999.0) {  /* valid entry for slope */
  79.       cs_2th += (double) ((1.0 - slpi * slpi) / (1.0 + slpi * slpi));
  80.       sn_2th += (double) (2.0 * slpi / (1.0 + slpi * slpi));
  81.     }
  82.   }
  83.   sn_2th /= (double) n_sgl;
  84.   cs_2th /= (double) n_sgl;
  85.   if (fabs (cs_2th) < CS_LOLIM)
  86.     cs_2th = SIGN (cs_2th) * CS_LOLIM;
  87.  
  88.   th = 0.5 * atan2 (sn_2th, cs_2th);
  89.  
  90.   if (fabs (sn_2th) < SN_LOLIM) {
  91.     printf ("\n...no preferred direction!\n");
  92.     th = 0.0;
  93.   }
  94. #ifdef DEBUG
  95.   printf ("...sn_2th = %6.4lf, cs_2th = %6.4lf, th = %8.4lf deg\n",
  96.           sn_2th, cs_2th, th);
  97. #endif
  98.  
  99.   av_slope = (float) tan (th);
  100.  
  101.  
  102. /* evaluate order parameter */
  103.   order_parm = (float) 0.0;
  104.   for (i = 0; i < n_segm; i++) {  /* deviations from av slope */
  105.     if ((slpi = *(dirns + i)) < 999.0) {
  106.       if (slpi * av_slope == -1.0)
  107.         tg_th = 512.0;          /* perp orientation */
  108.       else
  109.         tg_th = (slpi - av_slope) / (1.0 + slpi * av_slope);
  110.  
  111.       cs2_th = 1.0 / (1.0 + tg_th * tg_th);
  112.       if (DIM == THREE)
  113.         s = 0.5 * (3.0 * cs2_th - 1.0);
  114.       if (DIM == TWO)
  115.         s = (2.0 * cs2_th - 1.0);
  116.       order_parm += (float) s;
  117.     }
  118.   }
  119.   order_parm /= (float) n_sgl;
  120.  
  121.   return (order_parm);
  122. }
  123.  
  124.  
  125. /*
  126.  * eccentr()
  127.  *   DESCRIPTION:
  128.  *     eccentr evaluates aspect ratio of a given cluster (SGL) of segments
  129.  *     as the ratio of the maximum longitudinal and transverse dimensions,
  130.  *     the latter being identical to the maximum segment length in the SGL
  131.  *   ARGUMENTS:
  132.  *     segm(struct Segm *) pointer to segment structure (see sgl_stat.h)
  133.  *     list(struct linklist *) pointer to linked list structure (see lldef.h)
  134.  *     len(float *) length of SGL
  135.  *     width(float *) width of SGL
  136.  *   RETURN VALUE:
  137.  *     aspect ratio of SGL (length/width) 
  138.  */
  139.  
  140. float
  141. eccentr (struct Segm *segm, struct linklist *list, float *len, float *width)
  142. {
  143.   struct Segmtype *cSegm;
  144.   float del_x, del_y;
  145.   float d_min, d_max;
  146.   float csegm_len, max_segm_len;
  147.  
  148.  
  149. /* length */
  150.   llhead (list);
  151.   cSegm = (struct Segmtype *) list->clp->item;
  152.   d_min = cSegm->dij;
  153.   lltail (list);
  154.   cSegm = (struct Segmtype *) list->clp->item;
  155.   d_max = cSegm->dij;
  156.   *len = d_max - d_min;
  157.  
  158. /* width */
  159.   max_segm_len = (float) 0.0;
  160.   llhead (list);
  161.   do {
  162.     cSegm = (struct Segmtype *) list->clp->item;
  163.  
  164.     del_x = (float) ((segm + cSegm->segm_ind)->ptF.x - (segm + cSegm->segm_ind)->ptO.x);
  165.     del_y = (float) ((segm + cSegm->segm_ind)->ptF.y - (segm + cSegm->segm_ind)->ptO.y);
  166.     csegm_len = (float) sqrt (del_x * del_x + del_y * del_y);
  167.     if (csegm_len > max_segm_len)
  168.       max_segm_len = csegm_len;
  169.  
  170.   } while (llnext (list) == True);
  171.   *width = max_segm_len;
  172.  
  173.   return (*len / (*width));
  174. }
  175.  
  176. /*
  177.  * find_ex_hist()
  178.  *   DESCRIPTION:
  179.  *     find_ex_hist constructs a histogram of SGL (cluster) aspect ratios
  180.  *     (ex refers to eccentricity for which the aspect ratio is a crude measure)
  181.  *   ARGUMENTS:
  182.  *     n_segm(int) number of SGLs to calculate histogram for
  183.  *     asp_ratio(float *) array of aspect ratios for SGLs
  184.  *     ex_hist(int *) output histogram
  185.  *     bw(double) bin width of the histogram
  186.  *     ar_base(double) minimum aspect ratio (base)
  187.  *   RETURN VALUE:
  188.  *     none 
  189.  */
  190.  
  191. void
  192. find_ex_hist (int n_segm, float *asp_ratio, int *ex_hist, double bw, double ar_base)
  193. {
  194.   int i, ibin;
  195.  
  196.   for (i = 0; i < n_segm; i++) {
  197.     if (*(asp_ratio + i) != 0.0) {
  198.       ibin = 1;
  199.       while (*(asp_ratio + i) > ar_base + ibin * bw)
  200.         ibin++;
  201.       *(ex_hist + ibin - 1) += 1;
  202.     }
  203.   }
  204. }
  205.  
  206.  
  207. /*
  208.  * find_lsz_hist()
  209.  *   DESCRIPTION:
  210.  *     find_lsz_hist constructs a histogram of SGL (cluster) longitudinal sizes (no of segm)
  211.  *   ARGUMENTS:
  212.  *     n_segm(int) number of SGLs to calculate histogram for
  213.  *     lsize(int *) array of longitudinal for SGLs
  214.  *     lsz_hist(int *) calculated histogram for SGLs
  215.  *     nbins(int) currently not used
  216.  *     lsz_base(int) base of the histogram
  217.  *   RETURN VALUE:
  218.  *     none 
  219.  */
  220.  
  221. void
  222. find_lsz_hist (int n_segm, int *lsize, int *lsz_hist, int nbins, int lsz_base)
  223. {
  224.   int i, ibin;
  225.  
  226.   for (i = 0; i < n_segm; i++) {
  227.     if (*(lsize + i) != 0) {
  228.       ibin = 1;
  229.       while (*(lsize + i) >= lsz_base + ibin)
  230.         ibin++;
  231.       *(lsz_hist + ibin - 1) += 1;
  232.     }
  233.   }
  234. }
  235.  
  236.  
  237. /*
  238.  * zero_moment()
  239.  *   DESCRIPTION:
  240.  *     zero_moment evaluates zero order moment: x_vert and y_vert are assumed to represent
  241.  *     a closed polygon (last vertex assumed identical with first vertex)
  242.  *   ARGUMENTS:
  243.  *     nv(int) number of vertices in polygon
  244.  *     xv, yv(int *) (x,y) points of vertices
  245.  *   RETURN VALUE:
  246.  *     zero order moment of polygon 
  247.  */
  248.  
  249. int
  250. zero_moment (int nv, int *xv, int *yv)
  251. {
  252.   register int i;
  253.  
  254.   double xim1, xi, yim1, yi;
  255.   double m00, m00_sum;
  256.  
  257.   m00 = m00_sum = 0.0;
  258.   for (i = 1; i <= nv; i++) {
  259.     xim1 = (double) (*(xv + i - 1));
  260.     xi = (double) (*(xv + i));
  261.     yim1 = (double) (*(yv + i - 1));
  262.     yi = (double) (*(yv + i));
  263.  
  264.     m00_sum = 0.5 * (yi * xim1 - xi * yim1);
  265.     m00 += m00_sum;
  266.   }
  267.  
  268.   return ((int) fabs (m00));
  269. }
  270.  
  271.  
  272. /*
  273.  * find_area_hist()
  274.  *   DESCRIPTION:
  275.  *     find_area_hist constructs histogram of SGL (cluster) areas
  276.  *   ARGUMENTS:
  277.  *     n_segm(int) number of SGLs to calculate histogram for
  278.  *     area(int *) SGL array of areas
  279.  *     area_hist(int *) array of histograms filled
  280.  *     bw(int) bin width of the histogram
  281.  *     area_base(int) minimum aspect ratio (base)
  282.  *   RETURN VALUE:
  283.  *     zero order moment of polygon 
  284.  */
  285.  
  286. void
  287. find_area_hist (int n_segm, int *area, int *area_hist, int bw, int area_base)
  288. {
  289.   int i, ibin;
  290.   int up_bw;
  291.  
  292.  
  293.   up_bw = (int) bw + 1;
  294.   for (i = 0; i < n_segm; i++) {
  295.     if (*(area + i) != 0) {
  296.       ibin = 1;
  297.       while (*(area + i) > area_base + ibin * up_bw)
  298.         ibin++;
  299.       *(area_hist + ibin - 1) += 1;
  300.     }
  301.   }
  302. }
  303.